clearvars -except R bp bpSIM bgrid KAPPAS KAPPASIM BGpSIM BGgrid priceSIM priceSIM ySIM qtemp_big QSIM kappa_values PVD PVSIM PV_values yT_values G_big Def_func cte_debtfactor q NB NBG delta rbar omega ita EBGp Ebp EQ meank rhok sigmak NT BGp_big bp_cont q_cont ETau

 load Particle_endoPublic.mat
N=500000;
initialy=2008;
endy=2015; %Last Point in the Impact (sample goes to 2017)
first_point=1999;
window_push=initialy-first_point;
T=endy-initialy+1;
addpath('auxiliar')

%% Data
data_inputs= readtable('Output_PViip_NonPerf_9917.xls');
output=table2array(data_inputs(:,1));
pv_iip=table2array(data_inputs(:,2));
nonperf=table2array(data_inputs(:,3));
pub_iip=table2array(data_inputs(:,4));
spread_data=table2array(data_inputs(:,5));


nfa_data=pv_iip(window_push+1:window_push+T);
BGnfa_data=pub_iip(window_push+1:window_push+T);
y_data=exp(output(1+window_push:window_push+T));
pi_data=nonperf(1+window_push:window_push+T);
spread_input=spread_data(window_push+1:window_push+T);
% CENTERING
initialBG=pub_iip(window_push);
BGnfa_data=BGnfa_data-initialBG*ones(T,1);
nfa_data=nfa_data-initialBG*ones(T,1);

%% PArticle Filter

parfor BGi=1:NBG
    for i=1:NB
       % for inext=1:NBG
            b_VB(BGi,i,:,:,:)=reshape(squeeze(Ebp(BGi,i,:)),[6, 5,5]);
            BG_VB(BGi,i,:,:,:)=reshape(squeeze(EBGp(BGi,i,:)),[6, 5,5]);
            Q_VB(BGi,i,:,:,:)=reshape(squeeze(EQ(BGi,i,:)),[6, 5,5]);
            tau_VB(BGi,i,:,:,:)=reshape(squeeze(ETau(BGi,i,:)),[6, 5,5]);
       % end
    end
end    

bp_interpolate=griddedInterpolant({BGgrid,squeeze(bgrid(1,:)),yT_values,PV_values,kappa_values},b_VB,'linear','linear');
Q_interpolate=griddedInterpolant({BGgrid,squeeze(bgrid(1,:)),yT_values,PV_values,kappa_values},Q_VB,'linear','linear');
tau_interpolate=griddedInterpolant({BGgrid,squeeze(bgrid(1,:)),yT_values,PV_values,kappa_values},tau_VB,'linear','linear');
BGp_interpolate=griddedInterpolant({BGgrid,squeeze(bgrid(1,:)),yT_values,PV_values,kappa_values},BG_VB,'linear','linear');

%% Prepare to plot

maxB=max(b_SIM(:));
minB=min(b_SIM(:));
minBG=min(BG_SIM(:));
maxBG=max(BG_SIM(:));

bp_SIM=zeros(N,T);
bSIM=zeros(N,T);
BGp_SIM=zeros(N,T);
BG_SIM=zeros(N,T);
tauSIM=zeros(N,T);


QSIM=zeros(N,T);
spreadSIM=zeros(N,T);
bSIM(:,1)=b0;
BG_SIM(:,1)=BG0;

for i=1:N
for t=1:T
    bp_SIM(i,t)=min(max(bp_interpolate(BG_SIM(i,t),bSIM(i,t),ySIM(i,t),PVSIM(i,t),KAPPASIM(i,t)),minB),maxB);
      BGp_SIM(i,t)=min(max(BGp_interpolate(BG_SIM(i,t),bSIM(i,t),ySIM(i,t),PVSIM(i,t),KAPPASIM(i,t)),minBG),maxBG);
     QSIM(i,t)=min(max(Q_interpolate(BG_SIM(i,t),bSIM(i,t),ySIM(i,t),PVSIM(i,t),KAPPASIM(i,t)),1e-8),1);
       tauSIM(i,t)=max(min(tau_interpolate(BG_SIM(i,t),bSIM(i,t),ySIM(i,t),PVSIM(i,t),KAPPASIM(i,t)),1),0);
     spreadSIM(i,t)=(delta-delta*QSIM(i,t))/QSIM(i,t)-rbar;%rbar;%gSIM(i,t);%
    bSIM(i,t+1)=bp_SIM(i,t);
    BG_SIM(i,t+1)=BGp_SIM(i,t);
end
end




PubDebt=cte_debtfactor*BGp_SIM/mean_gdp;%./(Y_SIM_trans);
PVDebt=bp_SIM/mean_gdp;%./(Y_SIM_trans);

PubDebt_m=mean(PubDebt,1);
PVDebt_m=mean(PVDebt,1);
Q_m=mean(QSIM,1);
spread_m=(delta-delta.*Q_m)./Q_m-rbar;
nonperf_m=mean(PVSIM,1);
income_m=mean(ySIM,1);
kappa_m=mean(KAPPASIM,1);

time=[initialy:endy];



